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OPTIMALLY COUPLING THE KOLMOGOROV DIFFUSION, AND 
RELATED OPTIMAL CONTROL PROBLEMS 

KALVIS M. JANSONS AND PAUL D. METCALFE 
Abstract 

We discuss the optimal Markovian coupling before an exponen- 
tial time of the Kolmogorov diffusion, and a class of related stochas- 
tic control problems in which the aim is to hit the origin before an 
exponential time. We provide a scaling argument for the optimal 
control in the near field and use rational WKB approximation to ob- 
tain the optimal control in the far field, and compare these analytical 
results with numerical experiments. 

In some of these optimal control problems, in which the advec- 
tion velocity field is bounded, we show that the probability of suc- 
cess field agrees exactly with its leading-order asymptotic approx- 
imation in some areas of the plane, up to an undetermined multi- 
plicative constant. We conjecture a necessary and sufficient condi- 
tion for this behaviour, which is strongly supported by numerical 
experiments. 

1 . Introduction 

The coupling method is a probabilistic tool for studying the convergence to stationarity 
of a random process. By running two copies of the process at the same time, but imposing 
a dependence between the two copies to make them collide quickly, it is possible to obtain 
rigorous bounds on the convergence to stationarity of probability distributions (in the L\ 
norm). This technique — the coupling method f2l — is often used to prove convergence 
to stationarity by proving that collision is guaranteed in infinite time. 

This paper discusses the optimal Markovian coupling before an independent exponential 
time of two copies of the Kolmogorov diffusion. This study is motivated both by some 
recent work by 1 5 1 and by problems in fluid-dynamical mixing, though the latter application 
will be studied in depth elsewhere. 

Two copies of the Kolmogorov diffusion can be almost surely coupled in infinite time 
(see UJ for proof). To obtain a nontrivial optimization problem, it is necessary to add 
some degree of urgency to the problem. The simplest way of doing this is by requiring 
coupling before an independent exponential time; this time-limited problem was posed 
by |5J. The memory less nature of the exponential time means that the resulting coupling 
problem remains time-independent. This modified coupling problem can be handled as 
an optimal control problem; we use rational asymptotic methods, scaling arguments, and 
numerical experiments to derive optimal controls, and hence optimal coupling strategies. 

The fluid-dynamical motivation of the work comes from the simple model problem of 
mixing in a shear flow. The path of a passive tracer in a shear flow v{x) ~ yx is governed 
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by the Ito stochastic differential equations (SDE) 

dX,^dB,. 

(1) 

dY,^yX,dt + dW,, 

in which B and W are independent standard Brownian motions. For simplicity, we ne- 
glect diffusion in the streamwise direction, as is common in dispersion theory. Our pro- 
cess {X,Y) then becomes the Kolmogorov diffusion L6J, which is the simplest physically 
motivated model of thermally driven particle motion'. This does modify the mixing prob- 
lem near the origin, but the far-field behaviour of the mixing problem is qualitatively un- 
changed. We obtain a toy fluid-dynamical problem, but one which is of independent interest 

mni 

In ^we formulate our Markovian coupling problem as an optimal control problem, and 
study this optimal control problem numerically in |3l We see that the spatial dependence of 
the optimal control divides the plane into geometrically simple regions. We identify near- 
field and far-field regions of the plane and study the shapes of these geometrical features 
in these asymptotic limits in § ^4151 

Using scaling arguments, in ^we study the optimal control in the near field. By bal- 
ancing two modes of failure of the near-field motion, we obtain a surprising scaling law 
for the near-field behavior of the optimal control, which agrees well with the numerical 
results. In ^we use WKB analysis (i.e. the large deviation limit) to study our problem 
in the far field. We obtain a complete solution of the problem to WKB order, which again 
agrees well with our numerics. The near-field and far-field analysis, taken together, provide 
a complete qualitative picture of this optimal control problem. 

We discuss some variant versions of our basic optimal control problem in |6l which 
highlight noteworthy features of this class of control problems that can still be understood 
in the WKB limit. 

Finally, in ^ we give a complete description of the optimal coupling before an ex- 
ponential time of the Kolmogorov diffusion, and discuss the extension of these ideas to 
fluid-dynamical problems. 



2. Formulation 

The separation {X,Y) of two copies of the Kolmogorov diffusion is governed by the Ito 
SDE 

dX,=dB\'^-dBf\ 
dy, = YX,df, 

in which B^'^ and B^^' are both standard Brownian motions. (This also governs the separa- 
tion of two particles in a shear flow if the streamwise diffusion is ignored.) We are free to 
impose any dependence between the driving Brownian motions, and we do this in such a 
way as to drive the separation {X,Y) to the origin quickly (see HI). 

We now impose a Markovian dependence between the two driving Brownian motions 
b''^ and b'^'. By a convexity argument, the local optimization problem can be solved by 
requirmg dBj'' = ±dBp',and our problem reduces to a control problem in two dimensions. 

With little cost, we consider a richer class of control problems in which the velocity in 
the y direction, v(x), is an antisymmetric function of x (subject to regularity requirements), 

M.e. Brown 's Brownian motion 
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although it should be noted that only the linear case can be interpreted as an optimal cou- 
pling problem. We therefore choose, for each point {x,y), o{x,y) £ {0, \/2D} in such a 
way as to shepherd the process {X,Y) defined by the Ito SDE 

dX, = a{Xr,Y,)dBr, 
dY, = v{Xr)dt 

into the origin, where B is a standard Brownian motion. Our control here is very weak; we 
have freedom only to switch the cross-stream diffusion on or off, and it is not immediately 
obvious that we have sufficient control to steer the particle into the origin with non-zero 
probability. It is clear that almost all paths that eventually hit the origin must have infinite 
winding number about the origin. 

A simple — degenerate — coupling strategy allows us to bring the path within e of the 
origin, for arbitrary e > 0. Allow the particle to diffuse until it hits the line x = e in y < 0. 
The particle then hits y = at x = e. Using a sequence of such steps, the particle can then 
be made to hit the origin with probability 1 (see 1 1 1 for proof)^. 

To remove this degeneracy it is necessary to limit the time allowed to the particle. A 
simple and natural way to do this is to mark the path at an independent exponential time 
of rate X, and to require the particle to hit the origin before its path is marked. This is also 
in keeping with our model of mixing in a fluid-dynamical system; we require our paths to 
couple before they are swept apart into the bulk of the flow. 

Let Ho be the (possibly infinite) time that the particle first hits the origin, and r;^ be the 
time at which the particle path is marked. Now, we seek to calculate 

^(x,y) = supP^^')[//o<7i], (4) 

o 

where ' is the law of a particle started at {x,y), for a given control a, and = denotes a 
definition. Due to the antisymmetry of v, for all {x,y), <^{—x, ~y) — (^{x,y), and we use this 
symmetry freely; in particular all our numerical solutions are given in x ^ 0. 

To derive the optimal control, consider a path started at {x,y). We apply an arbitrary 
control a over the time interval [0,min(/!,//o)], and the optimal control O* over the time 
interval {mm{h,Ho),Ho]. Now the probability that this path hits the origin before being 
marked is 

P={l-Xh)E[^{Xh,Yh)]+o{h) 

^{l-Xh) (^^{x,y) + hv{x) ^ + 0^ +o{h) 



dy 2 dx^ 



[x,y)+h[v{x)^ + \(5^^-^]+o(h). 



We maximize this probabiUty by choosing 



a*(.,,)J^' (6) 
0, otherwise, 



^If one attempts to numerically compute a coupling strategy for the Kolmogorov diffusion that gives almost-sure 
coupling in infinite time, one obtains this degenerate strategy. E is set by the gridscale. 
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and we see that (|) satisfies the backward equation 



v(x) 



D 



8> 



= 0, 



(7) 



where = max(0,X). Note that we require the equality in the condition of (|6j to avoid 
pathologies, for example when there is an open set in which d(^/dx is constant. 



3. Numerical results 
We solve the backward equation Q by timestepping the partial differential equation 



8> 



(8) 



until a steady state is reached. We impose the initial condition 

^fn ^ (^'3') = (0,0) 

I 0, otherwise, 

and the boundary condition (|)(f,0,0) ~ 1. Imposing the symmetry <^{x,y) = i^{—x,—y), 
we only solve on jc ^ 0. Using a tanh mapping, we transform all occurrences of (— oo^oo) 
to (—1,1). For unbounded problems, we use operator splitting ||7l to avoid CFL-based 
stability restrictions: on each timestep we first solve the hyperbolic part of using linear 
interpolation, and second solve the diffusive part using forward Euler in time, with second- 
order finite differences in x. This is unnecessary for bounded problems, where we use 
forward Euler in time, first-order upwind finite differencing in y, and second-order finite 
differencing in x. 

The optimization selects a distinct region of the plane in which to turn the diffusion 
off, as seen in figure[2 The aim of this paper is to explain the distinctive shape of this no- 
diffusion region, which we sometimes refer to as region 2. Particles leave the no-diffusion 
region on the diffusion-on boundary y — a{x), and enter the no-diffusion region on the 
diffusion-off boundary y = b{x). In the far field, which we show in ^is governed by ray 
theory, we divide the diffusion region into two parts: regions 1 and 3. Ultimately successful 
particles started in region 1 move directly into the near field without crossing region 2. 
Ultimately successful particles started in region 3 cross region 2 in the far field before they 
reach the near field. We give a contour plot of the isf field for the case of figure^in figure|2 



3.1. The nature of successful paths 

A particle path that starts far from the origin, but hits the origin before being marked, 
has a distinctive shape which is observed in Monte-Carlo experiments. In its first wind 
about the origin the particle moves from the far field into the near-origin region. This 
motion is essentially ballistic, and is discussed further in ^ After arriving in the near- 
origin region the particle winds around the origin with large fluctuations, possibly making 
large excursions away from the origin. Eventually, the particle spirals into the origin in 
such a way that it never returns to the scale of a previous spiral, and each wind about the 
origin results in a huge decrease in the particle's distance from the origin. This motion is 
discussed in |3 The probability of this direct spiralling motion increases as the particle 
gets closer to the origin. 
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Figure 1: Boundaries of the portion of the no-diffusion region in x > for hnear flow 
v(x) = X, with X=l, and D = 2. The top boundary isy = a{x) and the bottom isy = b{x). 
The j agged lines for large x are due to both the decreasing density of our stretched grid and 
plotting artifacts. 




-3 I ^ ' ' ' ' 1 

0.5 1 1.5 2 2.5 3 

Figure 2: Contour plot of (|) in x ^ for linear flow v(x) = x, with X = \ and D = 2. The 
contours are at (|) = 3~", forn = 1 , 2, . . . , 7, from around the origin outwards. 
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Now, the probability of an unsuccessful wind, i.e. a wind that does not significantly re- 
duce the distance of the particle from the origin, must be a rapidly decaying function of the 
winding number, so that there is a non-zero probability of hitting the origin before mark- 
ing in the limit of infinite winding number In particular, this implies that the near-origin 
problem cannot be governed solely by Brownian scaling, as this would give a constant 
(limiting) probability of failure on each wind. This tension against pure Brownian scaling 
characterizes the dynamics of the near-origin region, and leads to the curious scaling laws 
described below. 

4. Near-field scaling law 

Near the origin, the position of the diffusion-off boundary is controlled by the balance 
between the two effects that lead to failure of the directly spiralling motion: either the 
particle path is marked before it hits the origin, or the particle fails to continue spiralling 
and returns to the scale of its previous wind. Because the relative decrease in scale on 
each wind becomes greater and greater, the estimate of the balance of these two sources of 
failure can be made (asymptotically in small x, or large winding number) by considering 
just one half-wind. 

Now consider one half- wind for a particle near the origin starting at a point {xi,yi) on 
the diffusion-off boundary. We suppose that the flow v{x) is antisymmetric, that v(x) > 
for X > 0, and that v(x) ~ |x|P for small x, where A ~ B means that A/B is (at worst) a 
slowly varying function^, which may be negative. In this case, the no-diffusion region lies 
in xy ^ 0. We assume that yi ~ |xi |", and we see that xiyi < 0. 

To obtain a consistent balance, we find that the dominant contribution to the probability 
of being marked on a half- wind comes from the motion in the no-diffusion region, and has 
probability Fi ~ 3;i/v(xi), for sufficiently small xi. For v(x) ~ |x|P, we find Fi ~ |xi 

Next, consider the particle motion from the diffusion-on boundary until it first hits the 
y axis. From Brownian scaling, the particle first hits the y axis at y2 — |xi Now the 
condition for the particle to continue to spiral into the origin is that it hits the diffusion-off 
boundary before it makes an excursion from the y axis with x displacement of scale xi. 
Thus this type of failure has probability F2 ~ x2/(|xi | + |x2|) ~ X2/X1, where X2 is given by 
y2 ~ |x2 r, with X23'2 < 0. So F2 ~ |xi | ' . The net probability of failure Fi + F2 in 

this half-wind is a minimum when Fi ~ F2, which implies that the scaling exponents are 
equal, giving: 

a2-a(p-l)-(p + 2) = 0. (10) 
The scaling law thus has exponent 

a = i(p-l + (p2 + 2p + 9)2). (11) 

For the Kolmogorov diffusion, we therefore predict a scaling exponent a = \/3. 

We attempted to determine the slowly varying function, or possibly constant, prefactor 
in the above scaling law, but were unable to adequately control the tails of some of the 
distributions we encountered (particularly in the y direction). 

Table n] compares the scaling (lll> with the results of the numerical solution of Q, in 
— 1 ^ X ^ 1 and —1 ^ y ^ 1, with killing boundary conditions. We compute the numerical 



^ A function f{x) is slowly varying if, for all Tl > 0, f(v\x) / f{x) ^ 1 as x [ 0. 
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Table 1: Scaling exponents in the near field for various power-law flows. 



p 




exact a 


numerical a 


1/2 


(a 


/4T-l)/4« 1.351 


1.372 


1 




\/3Ri 1.732 


1.760 


3/2 




/57 + l)/4Ri 2.137 


2.201 


2 




/T7 + l)/2Ri 2.562 


2.568 



exponents by a least-squares fit of 

log{-b{x))=c + mlogx (12) 

to the diffusion-off boundary b{x), where c and m are constants. 
This leads us to the following conjecture. 

Conjecture 1. Ifv{x) is antisymmetric, v(x) > Oforx> 0, and 

hml2^=P>0, (13) 
40 logx 

then the diffusion-off boundary b{x) satisfies 

lim i^lhH^)) = 1 (p _ 1 + (p2 ^ 2P + 9)^ ). (14) 
40 logx ^ 

Our numerical experiments do not exclude the possibility that if v(x) \x\^, then b{x) <x 
|x|", but this clearly cannot be shown definitively by numerical means. 

If, instead of the independent exponential time used elsewhere, we use a spatiaUy de- 
pendent marking rate ^(x) ~ |x|^ to limit the time, we can carry the same argument through 
to find 

a=i(|3-A-l + (p2 + 2p + 9+A2-2|3A + 2A)i/2), (15) 

although it is necessary to have p > A to obtain a non-degenerate solution. Numerical 
studies also support this scahng exponent, but this has not been verified to as high an 
accuracy as the above conjecture. 



5. WKB asymptotics in the far field 

As in the near field, the far-field scahng is governed by a balance of rare events. Here, 
we balance the probability that the particle path remains unmarked with the probabihty 
that the path has a sufficient diffusive motion to move to x = 0. If the particle takes a long 
time to move to x = 0, we increase the probability of the purely diffusive motion, but we 
also increase the probability that the particle path is marked, and vice versa. 

Throughout this section we assume that v is antisymmetric, continuously differentiable, 
and that v'(x) > for all x (where the ' denotes a derivative). This ensures that the no- 
diffusion region is connected, and its boundaries are functions of x. 

We derive the scaling of the most important regions of the far field by considering the 
X and y motion separately. First, we consider the motion in the x direction. The probabihty 
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density function for an unmarked Brownian particle in one dimension is 



p(x,xo,f) = (47tDf)" ' exp i^-^-^-Xtj , (16) 

and in the limit of large xq > 0, the dominant contribution of unmarked particles that hit 
X = is from paths of duration around fo, given by 

-^=Xto, (17) 
4Dfo 

balancing the probability of achieving the large diffusive stretch with the probability of 
remaining unmarked for a long time. (This implies a drift speed of V4Xd, which is exactly 
the drift of a Brownian motion conditioned on tending to infinity without being A--marked.) 
If we now consider the advection in the y direction with velocity scale V, we see that 
particles starting at {xo,yo), where 

yo ^ -Vto (18) 

will have both coordinates hit at roughly the same time. Thus particles starting in the 
neighbourhood of 

y~ -Vx/(A.D)'/2 (19) 

make the largest contribution to the flux of unmarked particles to the origin. The no- 
diffusion region is close to this region of the plane, so that it can control particles in this 
dynamically significant region. 

5.1. The meaning of 'far field ' 

These ideas can be refined to give the full leading-order behaviour of the diffusion-off 
and diffusion-on boundaries in the far field. We suppose that x scales as X, that the velocity 
scales as V, and thatj scales as XV / {4X0)^^^. RescaUng the backward equation Q gives 



8x2 



= 0, (20) 



where tildes denote nondimensional quantities and = {D /X)^I^X^^ . We now drop the 
tildes, and work with the non-dimensionalized equation j20> for the remainder of ^ Next, 
we seek a solution of the form (|)(x,3';e) = e.yjp{S{x,y\€) / and find that 5(x,y;e) sat- 
isfies the equation 



85y_ 2^ 

3x / 3x2 



i (21) 



on which we impose the boundary condition 5(0,0) — 0. We find the boundaries of the 
no-diffusion region y — a{x) and y — b{x) as asymptotic expansions in e, as e ^ 0: 

a{x) — aQ{x) + a\{x)z + a2{x)i^ + (9(8^), 

fc(x) = Z7o(x) + /7i(x)e + Z72(x)£^ + O(e^). 

The expansion in 8 is needed to split a double root. 

In the diffusion regions, we seek a solution of M\\ in WKB form: 

S{x,y-z)r^So{x,y)+Si{x,y)z + S2{x,y)i^ + --- . (23) 
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Terms with odd powers of e cannot be ruled out in region 3 due to e-scale features in region 
2. Substituting i23\ into (I2H and expanding, we find the equations 



^ ^ dy ^ dx dx ^ dx^ ^ \ dx ^ ^ 



at the first three orders. 



As is standard (see, for instance for more details), we solve the Hamilton-Jacobi 
equation (I24> by casting it as its equivalent set of Hamiltonian ordinary differential equa- 
tions. We first define a Hamiltonian 



H{x,y,px,py) = ~v{x)py + jp^ (27) 

from MA\ . and observe that the characteristics of the Hamilton-Jacobi equation MAI satisfy 
the Hamilton equations 



dx 


dH 


dx ~ 


dpx 


dy 


dH 


dx 


dpy 


dpx 


dH 


dx ~ 


dx 


Apy _ 


dH 


dx 





(28) 

V'{x)py, 



where x, y, px = dSo jdx and py = dSo /dy are considered as functions of the ray time x. The 
direction of the ray time has been chosen to ensure that information propagates away from 
the origin, and is in the opposite direction to the time of particle motion. Terms such as 
diffusion-off boundary and diffusion-on boundary refer to the underlying particle motion, 
and not to the motion of rays. We then find that the change in along a ray from {xs,ys) 
to {xe,ye) is 

A5o= ['^''^''\pxAx + pyAy), (29) 

but that S\ is constant on rays. 

For a ray that remains in the diffusion regions, we introduce the quantity 

r(xe,ye) 

5{xs,ys;xe,ye)= (pxdx + pydy), (30) 

Jixs.ys) 

where the integral is taken along the ray that runs from {xs,ys) to {xe,ye) with minimal 
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(31) 



ASq, so that that So{xi,yi) — S(0,0;xi,yi) in region 1. We note that 

OXe 

— {Xs,ys;Xe,ye) ^PyiXe), 

as 

dxs 

|^(x,.,y,;xe,ye) = -py{ls), 
ays 

where is and le are the ray time at the start and end of the ray, respectively. 

We extend the idea of rays to cover the no-diffusion region; here they have the Hamil- 
tonian —v{x)py = j, and we extend our definition of S with 



{xs,ys;xe,ye) = -Pxils), 



S{xs,ys;xe,ye) 



2v{xs) 



(32) 



where {xs,ys) and {xe,ye) are both in the no-diffusion region. Note that Xs = x^, on a ray in 
the no-diffusion region. However, S notation is not used for pieces of ray that cross from a 
diffusion region to a no-diffusion region or vice versa. 



5.2. Region 1 

In region 1 , we see that 

So{x,y) ^S(0,0;x,y), 
Sx{x,y)^Q, 

and we fit the boundary between regions 1 and 2 by requiring 



dx 



on this boundary. This condition gives 



dS 



,3^ 
3x2 



= 



^(x,ao(-'c)+) =0, 
ax 



(33) 



(34) 



(35) 



where the trailing + denotes a limit from above. This, not surprisingly, is also the condition 
for the action to be stationary under boundary perturbations. Using this condition, we see 
that 

px = {l-v{x)/v{x,)fl^ 



on a ray that hits the diffusion-on boundary at x = xo, and so 



a,{x) = -i v(^)(l-^ 



-1/2 



dt 



(36) 



(37) 



Now we can also compute the first-order correction to the diffusion-on boundary. On the 
diffusion-on boundary, and working to 0{e?'), the condition for the discriminant to change 
sign (I34> implies 



3^ 
3x3^ 



{x,aQ{x)+)ai{x) 



(x,flo(x)+) = 0. 



(38) 
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Next, we calculate the second-order derivatives in ( l38l . First, by differentiating (jSSj, we 
see that 

-^{x,aoix)+) + -^-^{x,ao{x)+)a'Q{x) =0. (39) 

Since the rays are vertical just above the diffusion-on boundary (px = 0, and py = 
— (2v(ji:))^'), we compute d^So/dxdy just above the diffusion-on boundary by differen- 
tiating along rays: 

|^(x,«o(x)+) = ^ (^) = -v'W/.>./v(x) = - , 
and so 

^(.,ao(x)+)=fl^(x)(^^j . (41) 

Now, we use the second derivatives ( 1401 14 1> in the perturbation discriminant J38> to find 
the first-order correction to the diffusion-on boundary: 

ai{x) = ^{x,aQ{x)+)(-^{x,ao{x)+)) = a'^ix) (^^j ■ (42) 

We choose the positive root of i42l to keep the discriminant positive in region 1 . 

Projecting our asymptotic expansion for a{x) into the class of decreasing functions, we 
obtain the approximation 

a{x) = ao{x) 1 - , (43) 



so that a{x) — a{x) — 0{e^). 
The action SQ{x,ao{x)) is then 



^aQ(x) 



5o(x,«oW) = - r^dt (44) 

We find, then, that for the power-law flow v{x) — sgn{x)\x\^, the first two terms in the 
asymptotic expansion of the diffusion-off boundary in x > are 

floW=-:tP+'B(il+r')/P 

, /o (45) 

a:W=xP+i/2[2(r^+r')B(il+r')]'/', 

and, just for interest, that 

5o(x,floW)=x(l+p-i)B(i,l+p-V2, (46) 
where B is the Beta function. 

5.3. Region 2 

To find the boundary b{x) of region 2 we find the points at which the discriminant J34> 
is zero. The action in region 2 is easily computed: 

S{x,y)^S{xMx)) + ^^^ 

\\_ (47) 

= 5o(x,aoW) + ^^^ff-^ + 0(e2). 
2v[x) 
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Substituting i4T\ into the discriminant, we find bo — ao, and then at first order we find 

-1 

(48) 

Here, we must choose the negative root to keep the discriminant negative in region 2; 
we see that bi = —ai. We also see that ( I48> has a root that gives the perturbation to the 
diffusion-on boundary; this is reassuring, as it shows that the analysis in region 1 is con- 
sistent with the analysis in region 2. Note that the symmetry between the diffusion-off and 
diffusion-on boundaries about y = ao{x) is not expected to extend to a2 and b2- 

5.4. Region 3 

The main purpose of ^is to determine the far-field approximations to the diffusion-off 
and diffusion-on boundaries, which can be done without reference to region 3. However, 
there are still some interesting features of the field in region 3, which are discussed here. 

By direct calculation, with error 0(£^), the action in region 2 is the analytic continuation 
of the action in region 1. Hence, the action in region 3 is, with error (9(e^), the analytic 
continuation of the action in region 1. Now the ray tubes in region 2 are perturbed (in 
position and cross-section) by 0{e^) due to the lack of diffusion, and they spend only 
(9(e) ray time, T, in region 2. This suggests that the true difference between the analytic 
continuation of the region 1 solution to region 3 and the true region 3 solution could be as 
small as 0{e^). 

The only remaining feature of region 3 is its boundary y — c{x) with region 1. This 
boundary is a caustic of the Hamilton-Jacobi problem (I24> . In region 1 , optimum rays head 
straight away from the origin, moving straight into positive x, whereas in region 3 optimum 
rays move into negative x before looping back into positive x. There is a transition region 
in the neighbourhood of the boundary of regions 1 and 3 where non-optimal rays exist; 
these rays are plotted for linear flow (v(x) = x) in figure|3] in which we launch rays from 
the origin to a point just inside the boundary of region 1. There are three rays from the 
origin to this point. Ray A first moves into x < 0, and has a large loop there before reaching 
X — 0. Ray B also moves into x < 0, but only has a small tight loop there. Finally, ray C 
moves directly into x > 0. Rays A and C give local minima of the action, but ray B gives a 
local maximum. By definition, the boundary y — c{x) between regions 1 and 3 is the curve 
on which the actions of rays of type A and type C are equal. Although figure|5]only shows 
the behaviour of rays in linear flow, we believe the qualitative behaviour of these rays is 
generic. 

For linear flow we explicitly determine the action in x ^ 0, and can therefore explicitly 
find this caustic. Define 

S±ix,y) = f(xV'C±(x,y)+T±(x,3;))±ix (49) 

where 

T±(x,3;) = (x2±63')'/2±x, (50) 

and then S{x,y) — 5_ (x,y) + O(e^) iny ^ c(x), and S{x,y) = S+{x,y)+0{e^) iny^ c(x). 
We find the caustic y = c(x) by solving S+{x,y) = S-{x,y) +0{e^), giving 

c(x) = -xV(4\/3) + 0(e2). (51) 



b]{x)=a'Q{x) 



2v(x) 
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Figure 3: Rays from the origin to a point just inside region 1, in linear flow (but ignoring 
the small effect from the no-diffusion region in x < 0, >■ > 0). The rays are plotted as solid 
lines, and the boundary between regions 1 and 3, and the line in region 1 on which non- 
optimal rays appear, are plotted as dotted lines. The caustic appears when the action along 
ray A and ray C are equal. Ray B is a local maximum of the action. 
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Figure 4: Comparison of numerical results and far-field asymptotics for the boundaries of 
the no-diffusion region, for linear flow v{x) = x, with A, = 1 , and D = 2. From top to bottom, 
the numerical diffusion-on and diffusion-off boundaries are plotted with solid lines, and the 
dimensional forms of the curves y = d{x), y = ao{x) and our asymptotic approximation for 
y — b{x) are plotted with dashed lines. The boundary, y — c{x), between regions 1 and 3 
is plotted with a dash-dotted line. The agreement between asymptotics and numerics is 
reasonable, given that e ~ 0.8, based onX =2. 

5.5. Far-field discussion 

Figure |4]plots the asymptotic results of ^for the far-field diffusion-off and diffusion- 
on boundaries, and gives a comparison with the numerically computed boundaries, for the 
case of a linear flow. Despite the large asymptotic parameter (e « 0.8), there is very good 
agreement between the asymptotic approximation of the diffusion-off boundary and the 
numerical boundary. The fit of the diffusion-on boundary is worse, although there is still 
reasonable agreement between the leading-order approximation ao{x) and the numerical 
boundary, which is between ao{x) and its first correction d{x). 

This large deviation calculation may equivalently be done in terms of a path integral for- 
mulation of the problem, in which the ray equations (I28> arise as the equations governing 
the classical path. 



6. 1 . Exact snap to far-field asymptotics 

Throughout this subsection, we assume that v is antisymmetric, continuous, v{x) > for 
X > 0, and without loss of generality that sup^v'(x) — 1. We further assume that there is a 
finite x,„ = inf{x : v(x) = 1 }. Under a local 'sharpness' condition of v in the neighbourhood 
of x,„ (given below), for sufficiently large \y\ the no-diffusion region degenerates to the 



6. Variant problems 



14 



Optimally coupling the Kolmogorov diffusion 








-0.5 



-1.5 



-2.5 



-3 







0.5 



1.5 



2 



2.5 



3 



Figure 5: A no-diffusion region degenerating to a line, for D — 2,X— 1, and with the flow 
field (|52|l. 

lines X = — sgn(3;)x,„, on which the field decays exponentially in \y\. This is seen in figure 
|5] which shows this phenomenon with the flow field 



Based on our WKB analysis, we conjecture a necessary and sufficient condition for 
the no-diffusion region to degenerate to lines x = —sgn{y)x,n in sufficiently large \y\. We 
simply ask whether a ray from the origin hits x — x„, with p^ — Q ata finite ray time (which 
is the condition for the ray to enter the no-diffusion region). This reduces to the condition 



for the no-diffusion region to degenerate to a line, and not otherwise. Based on numerical 
tests the above condition also appears valid for the case of bounded x, with reflecting 
boundary conditions, although the WKB argument given above is not directly applicable 
because of caustics due to boundary reflections. 

An example for which the above condition predicts no snap is v(x) = max(x(2 — x), 1) 
for X ^ 0, and this was numerically verified. This example is close to the marginal case. 

In some systems, and in some regions of the plane, this snap of the no-diffusion region 
leads to a curious behaviour in which the field (|)(x,y) agrees exactly with its leading-order 
far-field asymptotic solution, up to an undetermined multiplicative constant. We refer to 
this phenomenon as an exact snap. 

We now consider this exact snap in detail. From J20t . (^{x,y) ^ exp(— |e^^|y|) for all 
{x,y), based on the maximum advection velocity. Define y^ as the infimum of \y\ in the 




(52) 




(53) 
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degenerate part of the no-diffusion region. 
Define V|/ so that 



2 

then equation|20|implies for |v| > that 



^{x,y)=exp{-^e ^\y\)\\i{x,y), (54) 



vix)^ + ^£-'(1 +sgn(y)y(x))v = 0, (55) 



with discontinuities in d^]t/dx at x = ±x„, in respectively y < — y,„ and y > y,„. 

Under the above condition for a degenerate no-diffusion region, we find areas of the 
plane for which d\\t/dy = 0, namely y < —y,„ and x < x,„, and the corresponding area in 
y >0. Equation j55> then simplifies, and we are able to find the field in these areas of the 
plane, up to an undetermined constant ^{x,my,„). 

For y (x) = max(x, 1 ), for x ^ 0, we see that Xm — 1 , and that our WKB condition ( I53> is 
satisfied. In our exact snap areas, \55\ can be exactly solved in terms of Airy functions and 
exponentials. 

We test the hypothesis that d"^ /dy = in x ^ 1 for sufficiently large \y\ by computing 

For the data of figure El we find that ?v; (-oo, 1; -5.50, -7.46) = 5 x 10""^. (In figure 
|5] these dimensionless values of y correspond to dimensional values of —1.95 and —2.64 
respectively, based on X = 1 and V = 1. This gives e'* = 2.) It is difficult to get stronger 
numerical evidence for this in a doubly infinite system; our stretched grid coarsens away 
from the origin, and we choose y stations far enough from the ends of the grid to en- 
sure that the field is reasonably well resolved. For the values of y chosen here, the grid 
spacing in y is approximately 0.06. We find, from our numerical results, that the de- 
cay of the field along the filament of the no-diffusion region on x = — sgn(y)xm, \y\ > 
y,„ is exp(— 0.998 x jE^^y). Note, however, that the correlation coefficient between the 
grid points (x,) and (log(|)(x/, — 5.50)), in x ^ — 1, is 1 — 10^^. (We work in x ^ — 1 be- 
cause, for this problem, the predicted \\t field is particularly simple here.) The exponential 
(Ci exp(v^e^^x,)) (where Ci is a constant) is an excellent fit to the (|) field on x < — 1, 
y = —5.50. We also find that C2exp(V2e^^x+ jS-^^y) (where C2 is a constant) is a good 
fit to the (|) field over the area y < —ym, x < — 1. However, the decay at fixed y < —y^ of 
(|)(x,y) in X > 1 is (numerically) far from exp(— \/2e^^x), demonstrating that there is no 
exact snap in this area. 

A more convincing demonstration of this exact snap occurs in systems with bounded x 
and reflecting boundary conditions. In the corresponding system with x e [—1,1] and un- 
bounded y, higher resolution in y is possible, giving for example 9\C{—1, 1; —5.49, —6.62) 
< 10^' Here, the no-diffusion region degenerates ontox = — sgn(3') for sufficiently large 
\y\, and we are able to accurately resolve this boundary. 

The reason we believe that some systems have an exact snap is that this is the slowest 
decaying mode, so if the optimization can match to this solution it can do no better than 
this for larger \y\, given the amplitude of (|) at the matching y value. It appears that an exact 
snap does not occur in any area of the plane for which there exists a possible particle path 
to the X axis that does not pass through the no-diffusion region. 
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Figure 6; A gap in the no-diffusion region, for the flow field (|53, with D = 2 and A, = 1. 
The dashed fines show the ends of the interval (0.75,0.85) in which v is constant. 

6.2. Disconnected no-diffusion regions 

In the WKB calculation of ^ we assume that v is antisymmetric, continuously differen- 
tiable with v'{x) > 0, the last condition ensuring that the no-diffusion region is connected, 
significantly simplifying the calculation. If however we allow intervals \x\ G {^,^) where 
v'{x) = 0, these lead to gaps in the no-diffusion region over a range which includes (^,r|). 

An example for v with 

v'wJO' W ^ (0-75,0.85) ^3^^ 
I 1 , otherwise 

is given in figure |6l clearly showing a gap in the no-diffusion region due to the interval in 
which V is flat. Notice that the width of the gap in the no-diffusion region is larger than 
the width of the interval in which the velocity is flat, and this is a generic feature of such 
systems. Also generic is the filament of no-diffusion region that hangs down from the part 
of the no-diffusion region connected to the origin. This filament shortens as the gap size is 
reduced, and the y value at which the filament ends always roughly matches the y value at 
the bottom of the other side of the gap. 

If there are many such small gaps close together, the no-diffusion region is vertically 
striped, but its behaviour is almost the same as that of its envelope. This is because when 
a particle moves into the envelope, it either (i) enters a very x-restricted diffusion region 
which it passes through, (ii) enters a very jc-restricted diffusion region and later hits the 
adjacent no-diffusion region, or (iii) goes straight into the no-diffusion region. In all cases, 
the particle leaves the envelope with almost the same x. 

An extreme example of a striped envelope is when, apart from in the neighbourhood 
of the origin needed for success, the velocity is piecewise constant. In this case the no- 
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diffusion strips are of zero width. 



7. Discussion 

We are now in a position to clearly describe the behaviour of a particle that ultimately 
hits the origin before its path is A,-marked. For definiteness, we consider a particle in the 
unbounded shear flow v(x) — x, started at {Q,yo) with yo <0 and well into the far field. 
The particle first moves out into large jc > 0, to increase its flow speed towards the x axis. 
When the particle reaches y =yo/2 it turns around, and heads back towards x = 0. Near 
this maximum of x, the particle gains a small advantage by turning diffusion off, and hence 
truncating a quadratic maximum in the x motion, as seen in 13 The far-field motion (in the 
WKB limit) is little affected by the no-diffusion region, with corrections to the rays and the 
action first appearing at (9(e^) (and possibly higher), whereas the no-diffusion boundaries 
have corrections that scale as 8, caused by the splitting of a double root. The structure of 
the no-diffusion region for e = 1 is very similar to its structure in the far-field limit (i.e. 
small e). 

When the particle finally arrives near x = 0, it winds around the origin, with large fluc- 
tuations, until eventually it spirals directly into the origin. Controlling this rapid spiralling 
motion fixes the behaviour of the diffusion-off boundary near the origin, as seen in ^ In 
the near field the no-diffusion region is essential, as without it the probability of hitting 
the origin is zero. The near-field scaling laws are supported well by numerical experiments 
(including cases with power-law killing). A natural question is whether such scaling laws 
are possible for higher-order iterated diffusions, and whether it is possible to go beyond 
scahng laws to asymptotic expansions in the near field. 

The no-diffusion region is intrinsically governed by a balance between the probability 
that the particle path is marked, and the probability that a particle hits the origin at all. 
Considering the optimum no-diffusion region as a function of A,, we see that as A, increases, 
the probability of a particle starting at a fixed point hitting the origin marked or unmarked 
decreases. 

The natural fluid-dynamical extension of this work is to optimal coupling in a linear 
flow. The separation X between the two particles is governed by the Ito SDE 

dX,^A-X,dt + {I-R{X,))-dB,, (58) 

where B is Brownian motion, A is a traceless constant matrix, / is the identity matrix, and 
R is an orthogonal matrix. For each point, x, we choose R{x) in order to rapidly shepherd 
X into the origin, and hence achieve rapid coupling. We do this by forming the backward 
equation for the probability of successful coupling, analogous to Q. The local optimiza- 
tion problem, in this paper only a discrete choice, now becomes a continuum optimization 
problem, which may be solved by singular value decomposition fT|. The far field is stiU 
governed by ray theory, and similar techniques to those of ^ will still apply. However, 
the region corresponding to the no-diffusion region of this paper will have more structure, 
because the optimization is no longer a discrete choice. It may be possible to produce sys- 
tems of the form (I58> in which there is an intermediate scale resembling the near-field 
limit of ^ However, it is not clear if this effect will be important in real fluid-dynamical 
applications. 
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